/*
 #
 #  Files       : SquareMatrix.h
 #                ( C++ header file )
 #
 #  Description : The SmallMatrix Library
 #                ( http://code.google.com/p/smallmatrix )
 #
 #  Copyright   : Olivier Juan
 #                ( http://www.mas.ecp.fr/vision/Personnel/juan/ )
 #
 #  License     : CeCILL-C
 #                ( http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.html )
 #
 #  This software is governed by the CeCILL-C license under French law and
 #  abiding by the rules of distribution of free software.  You can  use,
 #  modify and or redistribute the software under the terms of the CeCILL-C
 #  license as circulated by CEA, CNRS and INRIA at the following URL
 #  "http://www.cecill.info".
 #
 #  As a counterpart to the access to the source code and  rights to copy,
 #  modify and redistribute granted by the license, users are provided only
 #  with a limited warranty  and the software's author,  the holder of the
 #  economic rights,  and the successive licensors  have only  limited
 #  liability.
 #
 #  In this respect, the user's attention is drawn to the risks associated
 #  with loading,  using,  modifying and/or developing or reproducing the
 #  software by the user in light of its specific status of free software,
 #  that may mean  that it is complicated to manipulate,  and  that  also
 #  therefore means  that it is reserved for developers  and  experienced
 #  professionals having in-depth computer knowledge. Users are therefore
 #  encouraged to load and test the software's suitability as regards their
 #  requirements in conditions enabling the security of their systems and/or
 #  data to be ensured and,  more generally, to use and operate it in the
 #  same conditions as regards security.
 #
 #  The fact that you are presently reading this means that you have had
 #  knowledge of the CeCILL-C license and that you accept its terms.
 #
*/

#ifndef _SMALLMATRIX_SQUAREMATRIX_H
#define _SMALLMATRIX_SQUAREMATRIX_H

#include <SmallMatrix/SmallMatrix.h>

NAMESPACE_BEGIN(SmallMatrix)

template <int N> class SquareMatrix: public Matrix_Base<SquareMatrix<N> >
{
private:
	typedef Matrix_Base<SquareMatrix<N> > Base;
public:
	template <int IdxR, int IdxC>
		reel& Access() {
			STATIC_ASSERT((IdxR<N)&&(IdxR>=0));
			STATIC_ASSERT((IdxC<N)&&(IdxC>=0));
			return Base::template dataAccess<IdxR*N+IdxC>();
		}
	template <int IdxR, int IdxC>
		const reel& Access() const {
			STATIC_ASSERT((IdxR<N)&&(IdxR>=0));
			STATIC_ASSERT((IdxC<N)&&(IdxC>=0));
 			return Base::template dataAccess<IdxR*N+IdxC>();
		}
public:
	INLINE SquareMatrix() : Base(){}
	INLINE SquareMatrix(const reel &v) : Base(v){}
	INLINE SquareMatrix(const reel tab[Base::SIZE]) : Base(tab){}
	INLINE SquareMatrix(const Vector<N>& left,const Vector<N>& right) {
		SquareMatrix<N>& ref=*this;
		const Matrix<1,N>& right_M = *reinterpret_cast<const Matrix<1,N>*>(&right);
		matrix_product_Unrolled<N>::template Operation<SquareMatrix<N>,Vector<N>,Matrix<1,N> >::compute(ref,left,right_M);
	}
	template <typename TT> INLINE SquareMatrix(const Matrix_Base<TT>& copy) : Base(copy){}
	INLINE SquareMatrix(const SquareMatrix<N>& copy) : Base(copy){}
	template <typename TT> INLINE SquareMatrix<N>& operator=(const Matrix_Base<TT>& copy) {
		Base& ref=*this; 
		return ref = copy;
	}
	INLINE SquareMatrix<N>& operator=(const SquareMatrix<N>& copy) {
		Base& ref=*this;
		return ref = copy;
	}
	INLINE reel& operator()(const int& idxm,const int& idxn);
	INLINE const reel& operator()(const int& idxm,const int& idxn)const;
	//template <typename T1,typename T2> INLINE SquareMatrix<N>& prod(const Matrix_Base<T1>& M1,const Matrix_Base<T2>& M2);
	//template <typename T1,typename T2> INLINE SquareMatrix<N>& add(const Matrix_Base<T1>& M1,const Matrix_Base<T2>& M2);
	//template <typename T1,typename T2> INLINE SquareMatrix<N>& subst(const Matrix_Base<T1>& M1,const Matrix_Base<T2>& M2);
	//template <typename T1> INLINE SquareMatrix<N>& add(const Matrix_Base<T1>& M1);
	//template <typename T1> INLINE SquareMatrix<N>& subst(const Matrix_Base<T1>& M1);
	//template <typename TT> INLINE SquareMatrix<N>& mult(const Matrix_Base<TT>& Two);
	//template <typename TT> INLINE SquareMatrix<N>& div(const Matrix_Base<TT>& Two);
	INLINE SquareMatrix<N> trans()const;
	INLINE reel trace()const;
	INLINE SquareMatrix<N> inv()const;
	INLINE SquareMatrix<N> inv(reel& deter)const;
	INLINE reel det()const;
	INLINE SquareMatrix<N>& addDiag(const reel& x);
};
//! Read/Write Accessor
template <int N> INLINE reel& SquareMatrix<N>::operator()(const int& idxm,const int& idxn)
{
	assert((idxm>=0)&&(idxm<N)&&(idxn>=0)&&(idxn<N));
	SquareMatrix<N>& ref = *this;
	return ref[idxm*N+idxn];
}
//! Read/Only Accessor
template <int N> INLINE const reel& SquareMatrix<N>::operator()(const int& idxm,const int& idxn) const
{
	assert((idxm>=0)&&(idxm<N)&&(idxn>=0)&&(idxn<N));
	const SquareMatrix<N>& ref = *this;
	return ref[idxm*N+idxn];
}
////! Matrix-Matrix Product in which the result is a SquareMatrix
//template <int N> template <typename T1,typename T2> INLINE SquareMatrix<N>& SquareMatrix<N>::prod(const Matrix_Base<T1>& M1,const Matrix_Base<T2>& M2) {
//	SquareMatrix<N>& ref = *this;
//	for (int i=0;i<N;i++)
//		for (int j=0;j<N;j++) {
//			ref(i,j)=0;
//			for (int k=0;k<T1::COL;k++)
//				ref(i,j)+=M1(i,k)*M2(k,j);
//		}
//	return ref;
//}
////! Matrix-Matrix Sum in which the result is a SquareMatrix
//template <int N> template <typename T1,typename T2> INLINE SquareMatrix<N>& SquareMatrix<N>::add(const Matrix_Base<T1>& M1,const Matrix_Base<T2>& M2) {
//	SquareMatrix<N>& ref=*this;
//	for (int i=0;i<N;i++)
//		for (int j=0;j<N;j++)
//			ref(i,j)=M1(i,j)+M2(i,j);
//	return ref;
//}
////! Matrix-Matrix Substract in which the result is a SquareMatrix
//template <int N> template <typename T1,typename T2> INLINE SquareMatrix<N>& SquareMatrix<N>::subst(const Matrix_Base<T1>& M1,const Matrix_Base<T2>& M2) {
//	SquareMatrix<N>& ref=*this;
//	for (int i=0;i<N;i++)
//		for (int j=0;j<N;j++)
//			ref(i,j)=M1(i,j)-M2(i,j);
//	return ref;
//}
////! Matrix-Matrix Sum in which the result is a SquareMatrix
//template <int N> template <typename T1> INLINE SquareMatrix<N>& SquareMatrix<N>::add(const Matrix_Base<T1>& M1) {
//	SquareMatrix<N>& ref=*this;
//	for (int i=0;i<N;i++)
//		for (int j=0;j<N;j++)
//			ref(i,j)+=M1(i,j);
//	return ref;
//}
////! Matrix-Matrix Substract in which the result is a SquareMatrix
//template <int N> template <typename T1> INLINE SquareMatrix<N>& SquareMatrix<N>::subst(const Matrix_Base<T1>& M1) {
//	SquareMatrix<N>& ref=*this;
//	for (int i=0;i<N;i++)
//		for (int j=0;j<N;j++)
//			ref(i,j)-=M1(i,j);
//	return ref;
//}
//template <int N> template <typename TT> INLINE SquareMatrix<N>& SquareMatrix<N>::mult(const Matrix_Base<TT>& Two) {
//	SquareMatrix<N>& ref=*this;
//	for (int i=0;i<N;i++)
//		for (int j=0;j<N;j++)
//			ref(i,j)*=Two(i,j);
//	return ref;
//}
//template <int N> template <typename TT> INLINE SquareMatrix<N>& SquareMatrix<N>::div(const Matrix_Base<TT>& Two) {
//	SquareMatrix<N>& ref=*this;
//	for (int i=0;i<N;i++)
//		for (int j=0;j<N;j++) {
//			assert(Two(i,j)!=0);
//			ref(i,j)/=Two(i,j);
//		}
//	return ref;
//}
//! SquareMatrix Transpose
template <int N> INLINE SquareMatrix<N> SquareMatrix<N>::trans()const
{
	const SquareMatrix<N>& ref = *this;
	SquareMatrix<N> res;
	matrix_trans_Unrolled<N>::compute(res,ref);
	return res;
}
//! SquareMatrix Trace
template <int N> INLINE reel SquareMatrix<N>::trace()const
{
	const SquareMatrix<N>& ref = *this;
	return trace_Unrolled<N>::compute(ref);
}

template <int N> INLINE	SquareMatrix<N> SquareMatrix<N>::inv() const {
	SquareMatrix<N> res;
	reel deter = det();
	assert(deter != 0);
	const SquareMatrix<N>& ref = *this;
	inverse_Unrolled<N>::compute(res,ref,deter);
	return res;
}
template <int N> INLINE	SquareMatrix<N> SquareMatrix<N>::inv(reel& deter) const {
	SquareMatrix<N> res;
	deter = det();
	assert(deter != 0);
	const SquareMatrix<N>& ref = *this;
	inverse_Unrolled<N>::compute(res,ref,deter);
	return res;
}
template <int N> INLINE reel SquareMatrix<N>::det() const {
	const SquareMatrix<N>& ref = *this;
	return det_Unrolled<IndexListEnd,N,N,N>::compute(ref);
}
//! Add \a x on the Diagonal of the SquareMatrix
template <int N> INLINE SquareMatrix<N>& SquareMatrix<N>::addDiag(const reel& x)
{
	SquareMatrix<N>& ref = *this;
	addDiag_Unrolled<N>::compute(ref,x);
	return ref;
}

NAMESPACE_END

#endif
